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ABSTRACT 

With a series of numerical simulations, we analyze the thermo-hydrodynamical evolution of circumstellar 
disks containing Jupiter-size protoplanets. In the framework of the two-dimensional approximation, we 
consider an energy equation that includes viscous heating and radiative effects in a simplified, yet consistent 



form. Multiple nested grids are used in order to study both global and local features around the planet. 
By means of different viscosity prescriptions, we investigate various temperature regimes. A planetary 
mass range from 0.1 to 1 Mj is examined. Computations show that gap formation is a general property 
which affects density, pressure, temperature, optical thickness, and radiated flux distributions. However, 
it remains a prominent feature only when the kinematic viscosity is on the order of 10 15 cm 2 s _1 or lower, 
though it becomes rather shallow for 0.1 Mj perturbers. Around accreting planets, a circumplanetary 
disk forms that has a surface density profile decaying exponentially with the distance and whose mass is 
5-6 orders of magnitudes smaller than Jupiter's mass. Circumplanetary disk temperature profiles decline 
roughly as the inverse of the distance from the planet, matching the values measured in the gap toward 
the border of the Roche lobe. Temperatures range from some 10 to ~ 1000 K. Moreover, circumplanetary 
disks are generally opaque, with optical thickness larger than 1 and aspect ratios around a few tenths. 
Non-accreting protoplanets provide quite different scenarios, with a clockwise, i.e., reversed flow rotation 
around low-mass bodies. Planetary accretion and migration rates depend on the viscosity regime, with 
discrepancies within an order of magnitude. Coorbital torques increase as viscosity increases. For high 
viscosities, Type I migration may extend to larger planetary masses. Estimates of growth and migration 
time scales inferred by these models are on the same orders of magnitude as those previously obtained with 
locally isothermal simulations, both in two and three dimensions. 

Subject headings: accretion, accretion disks — hydrodynamics — radiative transfer — methods: numerical — 
planetary systems: formation 

1. Introduction 
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Since the very first attempt of modeling a protoplanet 
interacting with its primitive nebula by Miki (1982), two 
decades have passed. Seventeen years had to go by for tack- 
ling this problem with the strength of new and more powerful 
computational tools (Bryden et al. 1999; Kley 1999; Lubow, 
Seibert, & Artymowicz 1999; Miyoshi et al. 1999). As a conse- 
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quence of the mounting interest in extrasolar planets, boosted 
by eight years of uninterrupted new detections 3 , the scientific 
community has made an extraordinary effort, over the past 
four years, to analyze disk-planet interactions by means of 
numerical calculations. 

Yet, regardless of the processes accounted for in the imple- 
mented numerical models, limitations and restrictions remain 
an issue to deal with. As a matter of fact, numerical calcu- 
lations still lag behind the present knowledge of the possible 
physical effects which occur in a circumstellar environment 
containing massive bodies. Thereupon, simplifying assump- 
tions are always demanded. Nonetheless, steps forward have 
been made in many directions to reduce their number. 

Three-dimensional (3D) disks containing Jupiter-size ob- 
jects have been investigated by Kley, D'Angelo, & Henning 
(2001) and, including also magnetic fields, by Papaloizou & 
Nelson (2003), Nelson & Papaloizou (2003), and Winters, Bal- 
bus, & Hawley (2003). D'Angelo, Henning, & Kley (2002, 
hereafter Paper I) carried out global two-dimensional (2D) 
high-resolution simulations of Jupiter-mass as well as Earth- 
mass planets, resolving also the flow structure within the 
planetary's Roche lobe by applying a resolution enhancement 
strategy. This study has been lately extended to the 3D space 
by D'Angelo, Kley, & Henning (2003, hereafter Paper II) 
and by Bate et al. (2003), who used instead a single grid with 
variable step size. Nelson & Benz (2003a,b) performed 2D 
computations, considering the effects due to disk self-gravity. 
Morohoshi & Tanaka (2003) have investigated the interac- 
tions between a planet and an optically thin disk with local 
simulations in the shearing sheet approximation, in order to 
determine the one-sided gravitational torque acting on the 
planet. 

With the present work we intend to relax the widely 
adopted local-isothermal approximation, i.e., that of treating 
a disk as a system having a fixed temperature distribution 
which depends neither on time nor on any other system vari- 
able. Yet, we will preserve the global viewpoint which we 
believe to be crucial because disk-planet interactions can be 
highly non-local. 

Many of the sophisticated accretion disk studies adopt 2D 
or 1+1D schemes due to the difficulty of dealing with all the 
necessary ingredients in a full 3D space. Until very recently, in 
fact, radiation (e.g., Dullemond et al. 2002), convection (e.g., 
Agol et al. 2001), magnetic fields (e.g., Matt et al. 2002) , 
mixtures of gas and dust (e.g., Suttner & Yorke 2001), and 
chemical evolution (e.g., Markwick et al. 2002) have been 
treated via two-dimensional models. Furthermore, numeri- 
cal simulations generally involve parameter studies, because 
many physical quantities are poorly known. Thereby, with a 
restricted number of dimensions, one is at least able to acquire 
a spectrum of possible approximate solutions. 

In order to overcome the local-isothermal hypothesis and 
to maintain a global approach at the same time, we intro- 
duce an energy transport in 2D (r-tp) models. We include 



Data archives on extra-solar planets, with the latest infor- 
mation, can be found at the Extrasolar Planet Encyclopedia 
(http://www.obspm.fr/planets) and at the California & Carnegie 
Planet Search (http://exoplanets.org/). 



all the supposedly major causes responsible for generation, 
transfer, and loss of energy in low-temperature circumstellar 
disk environments. We take advantage of the small aspect 
ratio of the disk, and assume that all the radiation trans- 
port is effective only in the vertical direction. Such approx- 
imation works reasonably well in accretion disks, away from 
the boundary and surface layers (see, e.g., Pringle 1981). 
Though such an approach relies on the "slimness" assump- 
tion, it is as well appropriate in the local environment around 
a protoplanet. In fact, around there, the disk becomes even 
thinner when the planet's gravitational action is accounted 
for. Hence, the amount of energy transported by radiation 
in the vertical direction still overwhelms that transported 
horizontally. Although this represents a simplistic kind of 
thermal description, it allows to investigate processes which 
have not been considered thus far, i.e., the joint thermo- 
hydrodynamical evolution of disk-planet interactions. Fur- 
thermore, since these simulations benefit of a nested-grid nu- 
merical technique (Paper I, Paper II), though relying on 
the global approach we are capable of achieving resolutions 
large enough to accurately describe the system within the 
planetary's Roche lobe. 

The outline of the paper is the following. In the next sec- 
tion we introduce the physical formulation of the problem, 
focusing on the equation for the thermal energy density and 
related issues. Section 3 concerns the numerical method uti- 
lized to solve the energy equation along with a series of tests. 
In § 4 we specify the adopted parameters and other various 
details. Sections 5 and 6 describe, respectively, global and 
local properties of models with different masses and viscosity 
regimes. In § 7 we report the results for mass accretion and 
orbital migration. Our conclusions are presented in § 8. 

2. Physical Formalism 

As mentioned before, we are going to tackle the problem of 
energy transfer in an accretion disk with an embedded planet. 
Since solving the full set of equations in three dimensions 
and simultaneously achieving sufficient resolution around the 
protoplanet is computationally difficult at the moment, we 
restrict to two-dimensional disk models. 

The reason for this choice is two-folded. On the one hand, 
with affordable computational times and sufficiently high res- 
olutions, it allows to investigate the coupling between hydro- 
dynamics and thermodynamics both in the disk and nearby 
the embedded body. On the other, it permits to adopt a 
strategic assumption in writing the energy equation, i.e., that 
the horizontal energy transport can be neglected compared to 
the vertical one. In other words, in a 2D disk one can assume 
that radiation transport is only effective along the direction 
perpendicular to the equatorial plane of the system. 

Hence, let us describe the disk material through the 
Navier-Stokes equations for the surface density S, the linear 
momentum (density) E u r , and the angular momentum (den- 
sity) r E Utp , where we choose as working coordinate frame a 
cylindrical one {O; r,ip, z} . The origin is located at the cen- 
ter of mass of the star and the planet while the disk lies in 
the plane z = 0. The hydrodynamics equations have been 
already stated, by using the same notations, in Paper I. 
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As in most of the computations performed so far, we avoid 
dealing with the Poisson equation and rather assume that 
the gravitational field is determined only by the star and the 
planet, neglecting the self-gravity of the disk material. This is 
allowed by the circumstance that the Toomre Q-parameter is 
well beyond 1 in all of our models. Nelson & Benz (2003a, b) 
found that some effects indeed arise when disk self-gravity is 
taken into account. However, we note that they consider disk 
masses ten times as large as those considered here. 

Thereby, we rely on a smoothed point-mass gravitational 
potential, which well reproduces more realistic protoplanetary 
gravitational potentials (Paper II), 



<1> 



GM, 



y/\r — r p | 2 + 8 2 



(1) 



where r„ and r p are the radius vectors indicating the po- 
sitions of the star and the planet, respectively. As done in 
Paper II, we set the parameter 8 as a constant length (see 
§ 4 for details). 

The thermal structure of the fluid has been usually fixed 
through a power law of the distance r. Thereupon, the set 
of equations has been closed by adding an equation of state 
which connects the local gas (two-dimensional) pressure to 
the surface density by means of a local isothermal sound speed 



= h 



GM, 



(2) 



Thus, the Mach number of the flow is considered constant 
throughout the system's evolution, constrained solely by the 
disk aspect ratio h — H/r (e.g., Lubow, Seibert, & Arty- 
mowicz 1999; Miyoshi et al. 1999; Nelson et al. 2000; Tanaka, 
Takeuchi, & Ward 2002; Masset 2002; Nelson & Papaloizou 
2003; Bate et al. 2003; Nelson & Benz 2003b). Alterna- 
tively, Kley (1999) and Tanigawa & Watanabe (2002) did 
some simulations using a polytropic equation of state of the 
type P — K TP . However, the evolution of the thermal prop- 
erties of the system have not been generally considered. 

In the present calculations we use an ideal equation of state 
which directly ties the gas pressure P to the thermal energy 
density (energy per unit area) E 



P = (l-1)E, 



(3) 



where the adiabatic index 7 is a constant. Furthermore, if we 
suppose that disk material behaves as a perfect gas, then the 
temperature is 

where p is the mean molecular weight, mu is the hydrogen 
mass, and k the Boltzmann constant. The adiabatic sound 
speed is given by 



kT 
prnn 



(5) 



Equation (3) implicitly assumes that radiation pressure 
frad is negligibly small compared to gas pressure. Such 
hypothesis is connected to the relatively low temperature 
regimes we deal with (in fact P ra d oc T 4 ) and its validity 
was checked afterward. In all of the models presented here, 
the ratio P r ad/P g as never exceeds 10~ 4 . 



2.1. Energy Equation 

Equations for energy transport differ according to the pro- 
cesses that have to be included for a correct description of the 
energy budget of fluid elements. In our case, we suppose that 
a parcel of disk fluid can gain or lose thermal energy only 
because of flow advection, compressional work, viscous dissi- 
pation, and dissipative effects due to radiation transport. In 
this sense, rather than strictly treating the transfer of radia- 
tion in the disk, we will account only for the cooling effects 
that radiation causes. Then, the energy equation takes the 
following form 



BE „ 



(Eu) = — P V ■ u + T — A. 



(6) 



In equation (6) we indicated with T the dissipation function 
and with A the radiated energy. In two dimensions, either of 
these quantities is an energy flux. Both functions are always 
positive, thereby the former acts as a heating source, whereas 
the latter as a cooling term. In this study we do not consider 
convection because we do not expect it to be a major energy 
source (D'Alessio et al. 1998). Furthermore, we neglect irra- 
diation from the central star which is presumably an effective 
heating mechanism only in the upper layers of circumstellar 
disks (see, e.g., D'Alessio et al. 1998), whose effects however 
are reduced by small disk scale heights. More importantly, 
as a consequence of the gap formation and the protoplanet's 
gravity, circumplanetary material is shielded from stellar ra- 
diation, as it will be clarified later. 

The function T can be directly computed from the compo- 
nents of the viscous stress tensor and S rv (Mihalas 
& Weibel Mihalas 1999; Collins, Heifer, & van Horn 1998) 



T 



2vT, 



[Si 



+ — — (V-w) 2 . 



(7) 



The divergence term V • u in the equation arises because a 
three-dimensional definition of stress tensor S is adopted. In 
fact, the reduction from three to two dimensions is merely 
achieved by setting the vertical component of the flow field 
u z and the derivative operator d/dz to zero. 

In a three-dimensional disk, the last term on the right- 
hand side of equation (6) is equal to V • F , where F is the 
frequency-integrated radiation flux 



F = _16^R T 3 VT 
3 Kp 



(8) 



In the above equation or is the Stefan-Boltzmann constant, k 
is a frequency-integrated opacity coefficient and p is the mass 
density. As mentioned before, we suppose that the amount 
of energy transported by radiation in the vertical direction is 
much larger than that transported horizontally, i.e., \F r \ and 
\F V \ <C \F Z \. The validity of this statement holds as long as 
the vertical extent of the disk remains very small compared 
to the disk characteristic size in the other directions. Thus, 
in a two-dimensional cylindrical disk we have 



/+00 /H 
V-Fdz~ / 
- OO J —I 



dF z 
dz 



dz. 



(9) 



Since the vertical disk structure is not meant to be accounted 
for, it is possible to equate the pressure scale height to the 
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photospheric scale height and assume that all the radiation is 
liberated at z — ±H. Thereupon, equation (9) becomes 



dF z 
dz 



dz = F(H) - F(-H) =2F(H). 



(10) 



In 2D disks it is useful to measure the emitted flux by means 
of an effective temperature: F(H) — anT 4 s , therefore 



A = 2<7 R T e 4 ff . 



(11) 



The factor 2 indicates that in a disk radiation escapes from 
both of its sides. 

A simple relation between the midplane temperature and 
the emergent radiation flux can be found by writing 



2 T— 

Jo dz 



dz 



= 2[F(H) 

8 OR 

3 KpH 



- F(0)] 
[T 4 (0) -T 4 (H)] 



(12) 



In the inner parts of a circumstellar disk the inequality 
T 4 (0) > T 4 (H) generally holds (e.g., Bell et al. 1997; 
D'Alessio et al. 1998), hence equation (12) yields 



A : 



2 a R T 4 
(3/4) r> 



(13) 



in which we introduced the optical thickness t = n p H = 
ttkE (we generally adopt a Rosseland mean opacity). We 
notice that with the previous definition the total disk optical 
thickness is 2 r. From now on, the quantity T refers to the 
disk midplane temperature. 

Equation (13) represents a fairly good approximation when 
the medium is very optically thick, i.e., r 3> 1. This is in- 
deed the case in those regions of unperturbed accretion disks 
which we simulate, because E is large enough. Yet, because 
of the action of massive bodies, in our case deep density gaps 
form where material is very diluted. In addition there are 
zones, where the disk spirals interact with the circumplane- 
tary disk spirals, which can have very low densities (see Pa- 
per I, Fig. 12 and 13). Such conditions make equation (13) 
not always applicable. Hubeny (1990) found a more rigorous 
relation between the effective and midplane accretion disk 
temperature, which represents a generalization of the gray 
model of classical stellar atmospheres in local thermodynamic 
equilibrium. Following Hubeny's theory, in a circumstellar 
disk we can write 



A = 2 cm T 4 



3r \/3 £h 
IT + ~4~ + 4t 



(14) 



According to equation (14), the emitted flux is inversely pro- 
portional to r in optically thick disk portions, whereas it be- 
comes proportional to r in the opposite limit. The quan- 
tity £h is equal to the ratio between the Rosseland and the 
Planck mean opacities. This quantity can be also roughly 
interpreted as the ratio of total extinction (i.e., absorption 
plus scattering coefficient) to the pure absorption coefficient: 
«cxt = ftabs + n BC . In these computations, we set e H = 1 be- 
cause we checked that we can neglect radiation scattering (see 
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Fig. 1. — Rosseland opacity coefficient k as function of the 
temperature for three different values of the mass density p 
(from Bell & Lin 1994). The legend quotes the density mag- 
nitude in gcm~ 3 . 



§ 3.1). Thus, no distinction is made between the extinction 
and the absorption coefficients. 

Equation (14) mimics the radiative losses of both optically 
thick and optically thin media, and therefore is very suitable 
to our purposes. In fact, it has been successfully applied to 
many studies on accretion disks (Popham et al. 1993; Godon 
1996; Collins, Heifer, & van Horn 1998; Hure et al. 2001). 

2.2. Opacity Table 

In order to calculate the disk optical semi-thickness r, we 
adopt the opacity formulas derived by Bell & Lin (1994) (see 
Fig. 1), as an improvement of those by Lin & Papaloizou 
(1985). Eight temperature regimes are identified according 
to the dominant processes active in each of them. Contri- 
butions from dust grains, molecules, atoms, and ions are ac- 
counted for. Since we simulate a distance range where disk 
material can be relatively cold, dust opacity is crucial for an 
accurate estimation of radiative losses (Lin 1981). Because of 
this, Bell's opacity includes grain absorption as tabulated by 
Alexander, Augason, & Johnson (1989). 

However, for the sake of comparison we ran some test cases 
(see § 3.1) with a new opacity coefficient developed by Se- 
menov and collaborators (2002, private communication) and 
based on the improved grain opacity tables provided by Hen- 
ning & Stognienko (1996). They aimed at coupling gas and 
dust opacities, focusing on the temperature range proper to 
protostellar disk environments. 

In either case, the midplane temperature T and mass den- 
sity p have to be provided in cgs units. In turn k = k(T, p), 
has units of square centimeters per gram. Within the in- 
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finitesimally thin disk approximation, dynamics variables are 
vertically integrated. Thus, the mass density p is not directly 
available. It must be retrieved from the surface density and 
the disk semi-thickness H. For this purpose, we assume that 



_E_ 
2H' 



(15) 



In the next section we explain how the disk scale height H is 
calculated in a fashion such to account also for the gravita- 
tional influence of the planet. This guarantees a more refined 
and consistent modeling of the system. 

2.3. Disk Scale Height 

Inside a circumstellar disk it is natural to assume that 
material is in hydrostatic equilibrium along the vertical di- 
rection (Pringle 1981; Frank, King, & Rainc 1992). Let g z 
be the gravitational acceleration in the z-direction within the 
disk, then the hydrostatic equilibrium reads 



g z dz. 



(16) 



The factor 2 is necessary because the midplane pressure sus- 
tains both sides of the disk. 

When the temperature distribution is not coupled to the 
rest of the dynamical variables, equation (16) is used to ob- 
tain the disk pressure from the density and the pressure scale 
height, assuming that g z = GM»/(r 2 + z 2 ) and omitting the 
planet gravitational attraction. 

In the present study, we compute the pressure directly 
from the thermodynamical processes occurring in the gas 
(eq. [3]), thus we can use equation (16) to estimate the disk 
scale height H in a consistent fashion. This quantity, in fact, 
is needed to obtain the mass density from the surface density, 
as explained above. Since the left-hand side of equation (16) 
is equal to C3/7, we can integrate the right-hand side in order 
to get an implicit function of H. Furthermore, the effects due 
to the gravitational field of the protoplanet can be included 
in this evaluation. If we set q — M*/M p and s = r — r p |, 
equation (16) generates the implicit function 



H 2 



2r ' 2 ^/{s/r) 2 + (H/r) 2 27 \G M, 



(17) 



Although there are no constraints on the ratio s/r, because 
with nested grid computations that ratio can be very small, 
the aspect disk ratio H/r is smaller than one, by working 
hypothesis, otherwise the 2D geometry would not be consis- 
tent. Therefore, the second term on the left-hand side can be 
expanded in a binomial series in H/r (up to the its second 
power), obtaining H as an explicit function 



H 2 = r 



7 \GM, 



+ q 



(18) 



We note that, in the limit ? -> or s > r, equation (18) 
yields y^y H = c s /Qk, as in regular accretion disks. Close 
to the planet, i.e., in the circumplanetary disk where s <JC r, 
the above relation reduces to yf^H = c s / \J G M p /s' A , which 
resembles the previous limit because the circumstellar Keple- 
rian angular velocity is replaced by the circumplanetary one. 



Rearranging the terms in equation (18), a more compact ex- 
pression can be written 



H 2 = 1 - 

7 



Cs 



l + q 



(19) 



Equation (19) has a singularity at s = 0. This arises from 
the corresponding singularity in the gravitational potential 
$, which is overcome by introducing the smoothing length 
(eq. [1]). In our computations, for consistency reasons, the 
distance \/s 2 + S 2 substitutes s in equation (19). 

2.4. Artificial Viscosity 

Non-linear effects in wave propagation inevitably lead to 
shock formation. This is indeed the case in disk-planet sim- 
ulations. In ideal fluids shocks are mathematical disconti- 
nuities. Therefore, in finite-differencing schemes, they must 
extend only over one or two grid points. In order to provide 
the correct jump conditions, ahead and behind a shock front, 
dissipative terms have to be present in the equations. This 
is usually done by introducing a non-linear viscous pressure 
otherwise known as artificial viscosity. Since the pressure 
gradient is no longer proportional to the density gradient, 
shocks stronger than those observed in local-isothermal com- 
putations (see, e.g., Paper I) may develop. Therefore, the 
physical viscosity might not be sufficient to provide the cor- 
rect jump conditions across shock fronts. 

The most rigorous treatment of shocks in a multidimen- 
sional space, with a generic metric tensor, requires the def- 
inition of an isotropic viscous stress tensor Q (Winkler & 
Norman 1986), which can be written as (Mihalas & Weibel 
Mihalas 1999; Stone & Norman 1992) 



Q = 



Vti - -(V ■ u) I 



(20) 



being I the unit tensor. Here we choose to discard the off- 
diagonal, i.e., shear tensor components because they may lead 
to artificial angular momentum transport (Stone & Norman 
1992). The coefficient of artificial viscosity is defined as 



MQ 



-£ 2 E min(V •«,()), 



(21) 



where £ represents the length over which the shock is smeared 
out. This is usually fixed to the maximum grid spacing. The 
coefficient /j,q is positive only for compression and zero for 
expansion, so the artificial viscosity is large in shocks and 
negligibly small elsewhere. 

Since the artificial viscous tensor acts as a pseudo-pressure, 
it affects both momentum and energy equations through the 
terms V • Q and Q • Vit (meant as a tensor product), respec- 
tively. The aforementioned equations are updated, by using 
the correct components Qij, as explained in Stone & Norman 
(1992). One severe side-effect is that Q can reduce the viscous 
time step by a factor which is proportional to [C 2 V • u] 1 . 

3. Energy Equation Solver 

The general numerical method employed to solve the hy- 
drodynamical equations on a hierarchy of nested grids, ap- 
plied to simulations of disk-planet interactions, has been ex- 
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plained in Paper I. In contrast to previous calculations, ther- 
mal energy density E now appears as an independent vari- 
able. Here we outline how the source terms (right-hand side 
of eq. [6]) are dealt with. In the framework of the nested-grid 
scheme, whenever required, E is interpolated from a finer to 
a coarser grid and vice versa, according to the procedures 
outlined in Paper I (see also Paper II). 

The energy equation (eq. [6]) is solved by means of a multi- 
step operator splitting method. The first step takes care of 
the energy advection and this is done in the same fashion as 
the advection of the other dynamical quantities, i.e., through 
the van Leer's algorithm. Viscous dissipation and radiative 
cooling are treated by means of a predictor-corrector scheme, 
which is second-order accurate in time, where heating and 
cooling terms are considered simultaneously for better accu- 
racy and stability. 

During the course of the first experiments, it was discov- 
ered that, due to very high pressures and large values of the 
flow divergence (in the wake of the circumplanetary spirals) 
the compressional work could lower the thermal energy by a 
large amount. Consequently, the predictor-corrector proce- 
dure applied to this equation term occasionally occurred to 
be unstable, producing negative energy values. Therefore we 
decided to take advantage of equation (3) and use an analytic 
solution for updating the energy. 

In the framework of the operator-splitting approach, in 
order to correct the energy because of the gas compres- 
sion/dilation, we have to solve numerically the equation 

dE 

¥ = - pv -- (22) 

With the aid of the state equation, this becomes 

^ = -(7-PV.u. (23) 

Equation (23) has an analytic solution, given by E = 
Eq exp [—(7 — 1) At V ■ u], where At is a time interval over 
which the divergence V • u has not appreciably changed. 
Hence, thermal energy density can be updated as follows: 



E n+1 = E„ cxp [-(7 - 1) At n V • «„]. 



(24) 



We notice that, expanding the exponential of equation (24) 
in a Taylor series and keeping all the terms up to the second 
order in At n , one finds 



E n +\ — E n 
At~ 



-P„ V • U n 



+ -(7-l)^A(„(V-«„) 2 , (25) 

where P„ = (7— 1) E n . It is easy to prove that equation (25) is 
equivalent to a predictor-corrector scheme. The above equa- 
tion also shows that the second-order correction always in- 
crease the thermal energy. 

Although the procedure in equation (24) is not so general 
as the predictor-corrector algorithm, it is more accurate and 
unconditionally stable. Furthermore, contrary to what may 
seem at first glace, from the computational viewpoint it is 
faster than the other because it requires less mathematical 
operations. 



3.1. Solver Test 

We have tested the equation energy solver from both the 
numerical and the physical point of view. Here we present 
some of these tests. 

The first test is intended to check whether the equation of 
energy furnishes physically consistent results by comparing 
the computed temperature to that derived from some analyt- 
ical solutions of equation (6). If we deal with a stationary Ke- 
plerian disk then the energy budget simplifies enormously. In 
fact, in a pure Keplerian flow both energy advection and com- 
pressional work are negligibly small. Therefore equation (6) 
reduces to (e.g., see Pringle 1981) 



T — A = 0. 



(26) 



From the same hypotheses, it follows that the dissipation 
function can be written in a simple form 



T = T.v 



AO, 



(27) 



In fact, in a Keplerian flow, the terms (du r /dr) 2 and (u r /r) 2 
are both much less than (r 80.K / 'dr) 2 . If we additionally as- 
sume that the disk is optically thick then the emitted flux can 
be approximated to 



A — Za-nJcs 



(3/8) t' 



(28) 



as implied by equation (14), when r >• 1. Recalling that 
r = isS, then, equating equation (27) and equation (28) 



yields 



128 o n 



(29) 



Now we assume that disk material has an opacity which we 
can be cast in the form 

K = K T 2 cm 2 g _1 . (30) 

with K = 2 x 1CT 6 (Morfill, Tscharnuter, & Voelk 1985). 
Placing equation (30) in equation (29), one finds 



J27_^ 
128 or 



(31) 



Such expression allows a direct comparison with the temper- 
ature distribution obtained from simulations, whose setup in- 
volves equation (27), equation (28), and the opacity law (30). 

In order to carry out a comparison of this kind, we simu- 
lated an unperturbed disk (i.e., without any embedded body), 
with borders at 1 and 20 AU, surrounding a Solar-mass star. 
The disk mass is Mo = 0.03 M* and the kinematic viscos- 
ity is v — 5 x 10 16 cm 2 s _1 . Both the initial surface density 
and temperature are constant: E(t = 0) = 197 gcm~ 2 and 
Tit = 0) = 352 K. Once the system has reached a station- 
ary state, the surface density is expected to decay as l/\fr, 
because v is constant (Lynden-Bell & Pringle 1974). Indeed, 
this is what one can observe in Figure 2 (left panel) , where the 
azimuthally averaged, computed surface density (£) (crosses) 
is fitted by the power-law 



E = 300 



5 AU 



gem" 



(32) 
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r [AU] r [AU] 

Fig. 2. — Test simulation of an unperturbed disk in which equations (27), (28), and (30) are implemented. Left: Crosses 
represent the computed E, averaged along the azimuth. The solid line indicates the power-law E = So y r p /r (see text). Right: 
The analytical formula for the midplane temperature (eq. [31]) is indicated in the plot by a solid line. Crosses represent the 
computed average temperature. 



On the other hand, equation (31) states that T ~ E f2 K ~ 
1/r 2 , or more precisely, T = 104 (5 AU/r) 2 K. Figure 2 (right 
panel) shows how the calculated averaged temperature (T) 
(crosses) fits to this analytic estimation. 

For the sake of completeness, we repeated the same type of 
test in which, instead of the Morfill et al.'s relation (eq. [30]), 
we chose the Kramer's law k — 6.6 x 10 22 pT~ 3 5 cm 2 g _1 , 
as in Frank, King, & Raine (1992), obtaining the same good 
agreement between analytical and numerical results. 

Although not reported here, we also tested the complete 
algorithm against convergence stability, i.e., we checked that 
it always converges to the same solution whatever the initial 
conditions are. Two simulations were set up in which the 
full form of equation (6) is solved along with equation (7), 
equation (14), and the Bell's opacity coefficient (§ 2.2). The 
disk radial domain is the same as before, but now the disk 
mass is smaller (Mo = 0.01 M») and so is the viscosity (y — 
10 16 cm 2 s _1 ). The initial E is imposed to be constant and 
equal to 66 gem -2 for both models. As initial temperature, 
in one model T(t = 0) = 14.5 K, while in the other T(t = 
0) = 580 K. The "cold" as well as the "hot" model soon 
evolve toward a stationary state. In each of the cases, the 
solutions exactly match, which indicates that initial values 
play no role in the equilibrium solution attained after the 
system has relaxed. The relaxation time somehow depends 
on the initial E, whereas scarcely any dependency on T at 
t — is observed. 

Finally, we examine some implications related to the opac- 
ity tables. A simulation was performed with the same param- 
eterization as the "hot" model mentioned above. But instead 
of choosing Bell's opacity, this model is based on the new 
tables by Semenov and collaborators (2002, private commu- 
nication). The resulting surface density distribution, at equi- 
librium, is hardly distinguishable from that obtained with 
the "hot" model. Furthermore, not much difference is mea- 



o 




. 1 I i i i . i i i i . i i i i . i i i \ 

5 10 15 20 

r [AU] 

Fig. 3. — Comparison between two protostellar disk models 
which differ only on the choice of the opacity tables. The 
solid thin line is produced by the model run with Bell's opac- 
ity formulas (Bell & Lin 1994), whereas the thick line comes 
from the model executed with the new opacity tables by Se- 
menov and collaborators (2002, private communication). The 
temperature is slightly larger in this second case with respect 
to the first one. Discrepancies generally stay within the 20%. 



sured between the stationary averaged temperature profiles, 
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as demonstrated in Figure 3. The result is not surprising, 
since already Liu & Meyer-Hofmeister (1997) tested the in- 
fluence of different opacity tables on the vertical structure of 
accretion disks and found that the change in disk structure, 
due to an improved opacity coefficient, is hardly perceivable 
when compared to the uncertainties connected with the gen- 
eral disk parameterization. 

Hence, throughout these simulations we use the opacity co- 
efficient by Bell & Lin (1994). It is worthy to note that this 
opacity table is well tested and has been extensively adopted 
in accretion disk studies (see, e.g., Godon 1996; Bell et al. 
1997; Klahr et al. 1999; Armitage, Clarke, & Tout 1999; Pa- 
paloizou & Terquem 1999; Nomura 2002). 

Equation (14) contains the quantity £h which can be in- 
terpreted as the ratio of extinction to the absorption coef- 
ficient. However, the physical conditions of the protostellar 
environment are such that it is allowed to neglect radiation 
scattering and write K cx t — K a t> s , and hence en ~ 1. In this 
phase, we checked such approximation by repeating the cal- 
culation addressed in Figure 3 and setting en = ^cxt /ft a bs , 
in equation (14). The resulting temperature profile is hardly 
distinguishable from the thick line in Figure 3. Therefore, the 
hypothesis en = 1 is quite correct in our applications. 

4. Model Parameters 

Model parameterization deserves more attention in these 
computations than it did in the simulations where energy 
transfer by radiation is not taken into account. The reason 
for this resides in the nature of the opacity coefficient k, that 
is always in cgs units. Hence, lengths and masses are to be 
fixed in physical units and consequently outcomes will not 
be scale-free, as they were in other circumstances (see, e.g., 
Paper I; Paper II). 

As in previous simulations, we model a circumstellar disk, 
orbiting a one-solar-mass star (M* = 1 M@), whose radial 
borders are r m i n = 2 AU and r max = 13 AU. The initial mass 
enclosed in this domain is M D = 4.8 x 10~ 3 M„, i.e., 0.01 M* 
within 20 AU (D'Alessio et al. 1998). Since in stationary Ke- 
plerian disks, the energy budget is regulated by the balance 
between viscous dissipation and radiative losses (eq. [26]), the 
magnitude of viscosity may play an important role. Hence, we 
will cover various viscosity regimes. The bulk of the simula- 
tions were run with a constant kinematic viscosity coefficient, 
by setting v = 1.0 x 10 15 , 5.0 x 10 15 , and 1.0 x 10 16 cn^s" 1 
(Bell et al. 1997). Besides physical viscosity, we consider an 
artificial viscosity, which is active only in compression regions, 
according to the tensor formulation explained in § 2.4. The 
length C is chosen to be equal to the maximum grid spacing, 
on each grid level (Stone & Norman 1992; Ziegler 1998). The 
mean molecular weight is /i = 2.39 (Morfill, Tscharnuter, & 
Voelk 1985) and the adiabatic index 7 = 1.4. 

Figure 4 shows characteristic averaged quantities of non- 
perturbed, stationary disks with the aforementioned charac- 
teristics. To distinguish among the different viscosity (and 
temperature) disk regimes, models will be referred to as "Hot" 
(H), "Warm" (W), and "Cold" (C). 

To further explore the parameter space, two models were 
executed applying the Shakura-Sunyaev viscosity prescription 



v = ac s H, with a = 10 2 and 10 3 (models A2 and A3). 
Hence, in such models v is a function of space and time. 

An embedded planet is placed on a fixed circular orbit at 
r p = 5.2 AU. The azimuthal position of the planet is kept 
fixed (<p p — ir) by the employment of a rotating grid. In order 
to obtain reliable outcomes in the two-dimensional geometry, 
we consider a planetary mass range extending from M p = 
33 Me to M p = 1 Mj or, in terms of mass ratios q = M p j M, , 
we have models with q = 1 x 10~ 4 ,2 x 10~ 4 ,5 x 10~ 4 , and 
9 = 1 x 10~ 3 . It has been recently demonstrated (Paper II) 
that, concerning gravitational torques and mass accretion, 3D 
models agree quite well with 2D ones in that mass range. 

However, since the disk thickness is sufficiently small, even 
for low protoplanetary masses (see Fig. 4), we go as far as to 
simulate an additional C- model with M p = 20 M©. Figure 5 
shows the disk semi-thickness in the planet's neighborhood 
for three selected masses. Indicating with Rn = r p (q/3) 1 ^ 3 
the planet's Hill radius, one can clearly see that the condi- 
tion H < Rn locally holds even in the smallest mass case. 
The figure also shows that protoplanets dwell in cavities and 
thus radiation from the central star cannot directly reach the 
surrounding matter. 

The smoothing length 5 is set to 2 x 10~ 2 Rn (see eq. [1]). 
All of the models are executed in an accreting and a non- 
accreting mode. When accretion onto the planet is allowed, 
the procedure outlined in Paper I is utilized. The accre- 
tion region has a radius K ac = 0.1 Rn, which should be short 
enough to let the whole accretion procedure be almost in- 
dependent of the removal characteristic time, as proved by 
Tanigawa & Watanabe (2002). 

Figure 6 illustrates the influence of the extent of the do- 
main's radial extent for Jupiter C-models. Simulations pro- 
vide quite consistent outcomes over the matching domain. 
Differences exist only close to the inner border, due to inflow 
boundary conditions. 

4.1. Initial and Boundary Conditions 

Initial conditions require the assignment of the three func- 
tions E(t = 0), T(t — 0), and u(t = 0) that we suppose to be 
axi-symmetric because we start from a stationary disk. 

The theoretical surface density profile of a stationary ac- 
cretion disk, with constant v, is a power-law of the radial 
distance r (Lynden-Bell & Pringle 1974), as also found nu- 
merically in § 3.1. Thus, as starting density distribution for 
our simulations, we set 

E(t = 0) = E M. (33) 

where E is determined by the disk mass value M D . 

Calculations carried out for different viscosity magnitudes 
imply different temperature regimes (see Fig. 4). According 
to them, the initial temperature distribution is prescribed. 
The behavior of the profiles in the top-right panel of Figure 4 
can be roughly reproduced by the power-law 

T(t = 0)=T {u) (34) 
The initial temperature at r — r p is given in Table 1 for three 
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Fig. 4. — Some characteristics of radiative, unperturbed and stationary disks, as function of the magnitude of kinematic 
viscosity (shown in the legend in cm 2 s _1 ). Radial boundaries are closed in order to achieve a strict stationarity by preventing 
mass losses. Top left: Surface density. Top right: Midplane temperature. Bottom left: Disk effective temperature (eq. [11]). 
Bottom right: Disk aspect ratio. 



values of v. Finally, the initial circumstellar flow u(t — 0) is 
a Keplerian one, corrected for the grid rotation. 

As for the boundary conditions, all models are run with a 
partially open inner radial border and a reflective outer one. 
Thus, matter is free to flow out of the computational domain, 
at r = r m i n , but the opposite is not allowed. This is the same 
expedient invoked in Paper I and Paper II to artificially 
mimic the mass accretion toward the central star and avoid 
spurious wave reflection at the inner domain edge, which is 
the closer to the planet. The flow field is Keplerian both at 
r min and r max : u = [0, r (Qk - fi P )]. 



4.2. Numerical Specifications 

The basic computational mesh is made of N r x N v = 
143 x 423 cells, while sub-grid patches are 64 x 64. All of the 
computations are based on a five-grid hierarchy. Thus, the 
finest resolution ranges from 1.3 x 10 -2 Rh (if M p = 1 Mj) 
to 2.9 x 1CT 2 Rh (if M p = 33 M®). Further details on nu- 
merical issues have been presented in Paper I and Paper II. 

Despite we have not implemented a sophisticated energy 
transport treatment, these simulations take between 35 and 
40% longer than those discussed in Paper I, for equal size 
grid hierarchies. 
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Fig. 5. — Because of the gravitational field of the embedded planet, the disk scale height may reduce considerably in its 
proximity, as prescribed by equation (19). This effect is shown here for models with different planetary masses: M p = 1 Mj with 
v = 1.0 x 10 16 cm 2 s- 1 (left); M p = 0.1 Mj with v = l.OxlO 16 cm 2 s" 1 (center); M p = 20 M© with v = 1.0 x 10 15 cm 2 s" 1 (right). 
For comparison, in the lowest mass case, Ru/r p = 0.027. 




Fig. 6. — Average surface density (left) and temperature (right) profiles belonging to two Jupiter-mass test models which differ 
only by the extension of the radial domain. The solid line refers to a model for which r m i n = 1.3 AU and r max = 20.8 AU 
(Md = 7.5 x 10~ 3 M»), whereas the dashed line represents a standard radial domain model [2, 13] AU (Md = 3.5 x 10~ 3 M„). 
The kinematic viscosity corresponds to that employed in C-models. 



5. Global Model Properties 

We now examine the behavior of global quantities for some 
selected models because the large scale look of disks with 
planets is, observationally, of extreme importance. Figure 7 
shows both the averaged surface density and temperature in 
Jupiter-mass computations with different kinematic viscosi- 
ties. We note that matter (left panel) is more uniformly 
distributed in the H- (solid line) and W-model (short-dash 
line), with respect to the C-model (long-dash line). In these 
cases, the density gap is not very deep and actually looks like 
a depression, because of the larger viscosity coefficient and 



the depletion of the inner disk. On the one hand, in fact, a 
high viscosity prevents the formation of wide and deep gaps. 
On the other hand, a large viscosity facilitates the inner-disk 
depletion, since Md oc i/E (Lynden-Bell & Pringle 1974). 
The disk depletion rate that we measure in these simulations 
ranges from few times 10 -8 to lO _7 M0yr _1 . 

In contrast, in the low viscosity regime (C-model) a well- 
defined and deep gap is carved in, where the density is more 
than an order of magnitude lower than it is in the other two 
models. The density peak visible in the middle of the gap 
is the signature of the circumplanetary disk. Since material 
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Table 1 

Initial Temperatures versus Viscosity and Disk Mass. 



v [cm 2 s x \ 


To[K] 


T max [K] 


T min [K] 


1.0 x 10 15 


50 


180 


13 


5.0 x 10 15 


110 


235 


22 


1.0 x 10 16 


150 


295 


26 



Note. — Values of the initial temperature as 
function the three values of the kinematic vis- 
cosity. Temperatures are sampled at r — r p — 
5.2 AU (T ), r = r max - 2.1 AU (T max ), and 
r = r min = 13 AU (T min ). Models with a- 
viscosity A2 and A3 are initiated as C-model. 




Fig. 7. — Average surface density (left) and temperature (right) in Jupiter-mass models after 240 orbital periods. By this 
evolutionary time, models have settled on a quasi-equilibrium state. The different profiles refer to different values of v, as 
indicated in the legends. The density gap, in the C-model (v — 10 15 cm 2 s _1 ), is more than an order of magnitude deeper 
(long-dash line) than it is in the other models. The temperature gap follows the same trend, though differences are smaller. 



removed from gap regions is disposed at short distances, due 
to the small damping length of the waves launched by the 
planet, gap shoulders have high densities. It has been argued 
that these density enhancements might trigger planet forma- 
tion. If so, low viscosity disks should be more efficient in 
doing that. 

Following the reduction of thermal energy in the low den- 
sity regions, a temperature gap accompanies the density gap 
(right panel, Fig. 7). In the C-model, temperatures drop be- 
low 20 K and it generally stays below 40 K. As expectable, 
the temperature gap has a counterpart in the pressure scale 
height, as seen in Figure 8. Such a deep trough represents 



an additional evidence that, as long as an inner disk exists, 
protoplanets are shaded against stellar radiation. 

Global views of both density and temperature 2D- 
distributions are displayed in Figure 9. Top panels illustrate 
the surface density, while the bottom ones show the tem- 
perature structure. From the density maps one can see that 
disk spirals get closer to each other toward the outer border 
of the domain. In fact, circumstellar spirals propagate at a 
velocity equal to the sound speed, which is proportional to 
VT (see eq. [5]). Thus, the propagation velocity reduces as 
r increases and therefore faster waves overtake slower ones. 
Past the gap region, temperature decays more rapidly in the 
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Fig. 8. — Average pressure scale-height (eq. [19]), normalized 
to the distance from the star, in Jupiter-mass models accord- 
ing to different viscosity values, as indicated in the legend (in 
cm 2 ^ 1 ). 



C-model (Fig. 7, right panel), hence this phenomenon is more 
evident. From the bottom panels one can also deduce that 
spiral perturbations in the thermal energy and mass densities 
are quite phase coherent and thus disturbances in the thermal 
structure are rather weak. 

Shifting toward smaller masses, the gap is gradually re- 
filled. This is visible in Figure 10, where C-models of different 
mass are compared. At M p = 0.1 Mj (i.e., M p = 33 M9), 
the gap is one order of magnitude as shallow as it is when 
M p —0.5 Mj. We will further address the issue of the gap 
structure in § 5.2, since gaps are considered promising features 
for detecting embedded planets. As regards the circumstellar 
disk spirals, only feeble traces are left, when M p — 0.1 Mj, 
both in the temperature and density distributions. The ab- 
sence of global disk features produced by an embedded body 
invalidates direct imaging observations as a detection tool 
in the intermediate-mass range, i.e., around a few times of 
Uranus' mass. 

Concerning the emitted energy, we can notice that the 
optical thickness of the models is generally larger than one, 
throughout the computational domain (Fig. 11, left panel). 
The averaged r stays above 1 also at the outer disk edge. 
This is in agreement with what was found for accretion disks 
around T Tauri stars (D'Alessio et al. 1998). The density and 
temperature gap reflects in a similar feature in the distribu- 
tion of r, as shown in Figure 11 (left panel). Nonetheless, 
even within the gap region, material is, on the average, op- 
tically thick (r > 1). The exception is represented by the 
C-model, where gap optical thickness reaches values around 
0.01 at ip = (p p + it (see Fig. 11, right panel). The most 



diluted portions of the gap, situated where the circumplane- 
tary disk merges into the gap, have values of r on the order 
of 10~ 2 -10~ 3 . 

Finally, we examine the emitted flux A, computed accord- 
ing to equation (14) (see Fig. 12). From the observational 
point of view, disk gaps could represent a probe for proto- 
planet detection. In fact, this is by far the most extensive 
imprint that a planet leaves on a circumstellar disks. Prospec- 
tive studies on observability of gaps due to Jupiter-mass pro- 
toplanets, have already been presented by Wolf et al. (2002); 
Steinacker & Henning (2003) , and speculations on a develop- 
ing gap around the T Tauri star TW Hya have been reported 
by Calvet et al. (2002). Here we do not intend to address the 
issue of whether gaps are really observable or to what ex- 
tent they are. However, we have to notice that a necessary 
but not sufficient condition is that they must be wide, deep, 
and with sharp edges in order for the flux emitted from this 
region to have the strongest contrast with respect to the sur- 
rounding environment (Rice et al. 2003). From the bottom 
panels of Figure 12, it appears clear that low kinematic vis- 
cosities would favor this kind of investigation. Disk spirals are 
probably too elusive to be captured by present-day ground- 
based instruments, as shown by the flux maps. We also ar- 
gue that it would seem rather unlikely to detect planetary 
masses smaller than Jupiter's, by means of these measure- 
ments. Right panels of Figure 12 display the emitted flux in 
case of a M p =0.1 JV/j model. The quantity A furnished by 
high-viscosity model (top) looks quite smooth and only some 
inhomogeneities appear in C-model (bottom). 

5.1. Models with a-viscosity 

So far we have considered disks with a constant v. How- 
ever, we have also investigated the effects of a different vis- 
cosity prescription, namely an a-type viscosity (Shakura & 
Sunyaev 1973). Two Jupiter- mass accreting models were 
run in this mode, one with a = 10~ 2 (model A2) and the 
other with 10 (model A3). Such values are representative 
of those recently estimated by Papaloizou & Nelson (2003) 
and Winters, Balbus, & Hawley (2003) who performed MHD 
disk simulations, with embedded Jupiter-size bodies, and 
could therefore evaluate the stress parameter a in a self- 
consistent fashion. Indeed, referring to the constant viscos- 
ity models discussed above, by making the azimuthal mean 
2-k/ (a) = / c s H dtp, we get a value of {a} ranging from 

3 x 10~ 2 in the inner disk to 4 x 10~ 3 in the outer disk. 
Specifically, C-models yield an overall trend of (a) which is 
the closest to 10~ 2 . And in fact outcomes from a-viscosity 
model A2 are not much different from those obtained from 
the Jupiter-mass C-model. The quantity (E) is very similar 
to that illustrated in left panel of Figure 7 (long-dashed line) 
whereas (T) is only slightly different (near the inner border) 
from the one in the right panel of the same Figure. 

Some more important deviations from what has been seen 
thus far are instead encountered in the A3-model, as displayed 
in Figure 13. The averaged surface density shows a significant 
hump (top-left panel, solid line), inside the gap region, not 
observed in the other viscosity regimes. This is not caused by 
the mass accumulation around the planet but it is produced 
by material lingering around Lagrangian points L4 and L5, 
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Fig. 9. — Two-dimensional surface density (top) and temperature (bottom) distributions for Jupiter-mass models. In the top 
panels, the color scale is linear; in the bottom ones the colors scale logarithmically. In this plot, E = 1CF 4 corresponds to 
32.9 gem -2 , whereas T — 1CP 2 equals 198 K. Left panels correspond to H-models whereas right panels refer to C models. 
Density distributions show clear evidences that the a well-defined gap needs a cold (i.e., low viscosity) environment to be 
established. The same consideration applies to the temperature distribution. 
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Fig. 10. — Left. Averaged surface density and temperature in C-models (y = 10 15 cm 2 s _1 ) with different masses: M p — 0.5 Mj 
(solid line), M p — 0.2 Mj (short-dash line), and M p = 0.1 Mj (long-dash line). Right. 2D-distribution of the surface density 
(top) and temperature (bottom) of a C-model with M p = 0.1 Mj. The density is scaled linearly, while temperature is scaled 
logarithmically. Conversion factors are: E = 10~ 4 — > 32.9 gem -2 , T = 10~ 2 — » 198 K. 
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Fig. 11. — Circumstellar disks are usually optically thick. This plot shows that this is generally the case even when a disk 
contains a Jupiter-mass body. The line sequence is the same as in Figure 7 (data refer to the same models). Left. Optical 
thickness averaged around the star. Only in the C-model, (t) drops below 1. Right. Profiles at (p ~ 0, i.e., 180° away from the 
planet location. The medium becomes very optically thin (r <C 1) only in the gap region of the C-model. 




Fig. 12. — Flux radiated (i.e., A, see § 2.1 for details) by H-models (top) and C-models (bottom). Left panels regard simulations 
of Jupiter-mass objects, whereas right panels describe models with M p — 0.1 Mj (M p = 33 M©). The color scale in both panels 
is logarithmic. In plot units, 1CP 5 means an energy flux of 4.7 x 10 4 ergcm _2 s _1 . For comparison, the surface flux of the Sun 
is 6.3 x 10 10 ergon^s -1 . 
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Fig. 13. — Le/t. Averaged surface density and temperature in the Jupiter-mass a-viscosity models A3 (a — 10 , solid line) 
and A2 (a = 1CP 2 , long-dash line). Right. 2D-distribution of the surface density (top) and temperature (bottom) from models 
A3. The density scales linearly and temperature scales logarithmically. Conversion factors are: E = 1CT 4 -> 32.9 gem , 
T = 1CT 2 -> 198 K. 
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as it appears from the ring-like traces in the middle of the 
density gap in the top-right panel Figure 13. This also appear 
to be a persistent feature, in fact no change is measured over 
the last 50 orbital periods. L4 and L5 points are equilibrium 
locations for the restricted three-body problem. Hence, one 
can expect that the thermal energy is larger near such points 
than anywhere else in the gap. In turn, viscous torques, which 
are enhanced by the locally increased kinematic viscosity are 
more efficient in contrasting gravitational torques. Thereby, 
material is more likely forced to evolve on the much longer 
viscous time scale t v i B ~ r p /u = r p /(ac B h). 

As for the temperature (Fig. 13, bottom panels), the 
largest discrepancies between models A2 and A3 occur to- 
ward the inner border. 

5.2. Gaps as Planet Signatures 

We have seen in the previous two sections that the gap 
structure drastically depends on the viscosity regime. The 
standard criterion to open a gap in a disk requires two con- 
ditions (Lin & Papaloizou 1985, 1993). The first is a thermal 
condition given by: 



Mp 

M* 



> 3 



-) -3 ft 3 , 
r 



(35) 



i.e., the Roche lobe must be larger than the disk scale height. 
Because of planet's gravity (see eq. [19]) and the low temper- 
atures established in gap regions, this condition is always ful- 
filled when M p > 0.1 Mj (see Fig. 8). The second condition 
concerns the requirement that tidal torques exceed viscous 
torques and reads 



M* 



(36) 



where 1Z = r 2 0,-k/v is the Reynolds number. The right hand 
side of equation (36) is proportional to u and in our case 
goes from Ri4x 10~ 3 (H-models) to w 4 x 10~ 4 (C- models). 
Thereby, we should not obtain gaps in any of the H-models. 
Figure 7 (solid line) decently agrees with this prediction, since 
only a trough is dug in H-model. However, condition (36) is 
partly violated in the low- viscosity regime (Fig. 10), because 
a clear gap is visible for the M p — 0.2 Mj case and a trough 
appears for M p = 0.1 Mj. In Table 2 we report the gap 
occurrence and depth for some of the models. 

One reason why the above criterion could not properly 
apply to our calculations lies in the fact that it was derived 
for simple polytropic disks, while we simulate also the thermal 
evolution of the system. As implied by the profiles in the top- 
left panel of Figure 13, it is even harder to predict the gap 
presence when the kinematic viscosity depends on the local 
fluid temperature. 

6. Circumplanetary Disks 

The first distinction about the circumplanetary flow has to 
be made according to whether the planet is accreting or not. 
Therefore, we shall carry out separate discussions, based on 
H- and C-models, as well as a-viscosity models A2 and A3. 



6.1. Accreting Models 

Figure 14 indicates that, inside of the Roche lobe, the gen- 
eral aspect of the flow circulation around 1 Mj and 0.1 Mj 
protoplanets resembles that obtained with local isothermal 
models in two dimensions (Paper I). However, specific char- 
acteristics of the flow do differ. A circumplanetary disk, ex- 
tending over a fair fraction of Roche lobe, can be identified 
for both masses and viscosity regimes. The mass of such 
disks around one Jupiter-mass protoplanets, within the 80% 
of the Hill radius, is 1.03 x 10~ 5 Mj in the H-model and 
5.36 x 10~ 6 Mj in the C-model. As regards calculations with 
M p = 0.1 Mj planets, the mass measured in the circum- 
planetary disk (or sub-disk, for brevity) of the C-model is 
4.03 x 10~ 6 Mj, while the amount is nearly two times less in 
the H-model. While the flow dynamics in the Jupiter-mass 
model with a — 10 -2 (left panel in Fig. 15) is similar that of 
the C-model, the model with a = 10~ 3 presents some pecu- 
liarities. As shown in right panel of Figure 15, the size of the 
circumplanetary disk is smaller than the Roche lobe. More- 
over, the streams of matter associated with the spiral waves 
are narrower. The sub-disk mass measured in these cases is 
3.9 x 10~ 6 Mj (A2) and 1.5 x 10~ 6 Mj (A3). 

The azimuthal average of the surface density distribution 
around the planet, inside the Roche lobe, can be fitted by a 
relation of the type (see Fig. 16): 



(E) ~ (E)o exp 



(37) 



where s/Rh < 1 and larger than a threshold length (either 
0.1 or 0.2 Rh). The values of (E)o and oi are reported in 
Table 3. 

The presence of spiral features in the density distribution 
is an indication that the circumplanetary flow is Keplerian- 
like. Indeed, decomposing the velocity field u into the in-fall 
velocity (toward the planet) Wi n and the rotational velocity 
(counter-clockwise around the planet) u> rot , it turns out that 
the former is more than an order of magnitude less than the 
latter. In Jupiter-mass models, the rotational velocity within 
0.5 Rh of the planet drops approximately as s -0 ' 6 , nearly in- 
dependently of the viscosity regime investigated in these com- 
putations. For 0.1 Mj models, this ratio declines at a slightly 
higher rate. Consequently, for s £ [0.1,0.5] Rh, the ratio of 
TOrot to the local Keplerian velocity is always larger than 0.7. 

Comparing the right panels (C-models) of Figure 14 to 
Figure 13 (top and middle right panels) in Paper I, one can 
clearly see that spiral perturbations are less intense, as we 
show below, and more open. The latter circumstance is re- 
lated to the lower values of the Mach number in the circum- 
planetary flow, which governs the inclination of the spiral 
wave with respect to the direction of the rotational motion. 
In fact, in local isothermal models A^iso — \J q (r p /s)/h. If 
M p = 1 Mj, Miso drops from ^ 8, at s — 0.1 Rh, to 2, at 
s — IRh- As comparison, in the circumplanetary disk dis- 
played in the upper-left panel of Figure 14, M lies between 
2.4 and 1.4, whereas the low viscosity Jupiter-model (upper- 
right panel) provides values between 3 and 1.5. Models A2 
and A3 yield similar values. The main reason for this differ- 
ence resides in the larger value of the sound speed. Hence, 
perturbations can travel faster and are less distorted by the 
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Table 2 
Gap Occurrence and Depth. 



Mp/Mj 




C-MODELS 






H-MODELS 




Acc. 


Non-Acc. 


N.F. 


Acc. 


Non-Acc. 


N.F. 


1.0 


< 0.01 


< 0.01 




0.6 


0.8 


X 


0.5 


0.06 


0.2 




0.9 


NG 


X 


0.2 


0.3 


0.4 


X 


NG 


NG 


X 


0.1 


0.6 


NG 


X 


NG 


NG 


X 


0.06 


0.8 


NG 


X 






X 



Note. — Ratio of the minimum density in the gap to the density at 
the gap's taller shoulder (rounded to the first significant digit). Values 
are recovered from the azimuthal average of S around the star. The 
letters "NG" stand for "No Gap" . The symbol "X" appear whenever 
at least one of criterion (35) and (36) is not fulfilled (N.F.). 



Table 3 

Fit Parameters for the Averaged Surface Density. 



M p /Mj - 




C-MODELS 


H-MODELS 


(E) [gem 


-'] ai Range 


(S)o [gem -2 ] ai Range 


1.0 
0.5 
0.2 
0.1 
0.06 


2.92 x 10 2 
1.59 x 10 2 
2.32 x 10 2 
1.99 x 10 2 
1.47 x 10 2 


-4.4 [0.2,1.0]i? H 
-2.0 [0.2,1.0]i? H 
-1.5 [0.1, 1.0] ifa 
-0.9 [0.1,1.0]i? H 
-0.4 [0.2, 1.0] .Rh 


2.14 x 10 2 -2.3 [0.2, 1.0] -R H 
1.35 x 10 2 -1.3 [0.2,1.0]i? H 
8.57 x 10 1 -0.5 [0.1, 1.0] -R H 
6.75 x 10 1 -0.2 [0.2, 1.0] -R H 






A3-MODEL 


A2-MODEL 


1.0 


1.80 x 10 2 


-7.5 [0.2,0.8]i? H 


2.36 x 10 2 -5.0 [0.2, 1.0] .R H 


Note.— 


-Parameters that enter equation (37). 


This represents a linear best-fit of the 



logarithm of the averaged density (£) over the specified interval. 
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Fig. 14. — Density distribution and velocity field around M p = 1 Mj {top) and M p = 0.1 Mj (bottom) models. Right panels 
refer to H-models (v — 10 16 cm 2 s _1 ) whereas left ones refer to C-models (y — 10 15 cm 2 s _1 ). In each panel, color scales are 
linear. In the units used in the color bar, 1CP 4 corresponds to £ = 32.9 gcm~ 2 . 



background motion of the flow. In the Jupiter-mass cases, 
illustrated in Figure 14, the azimuthally averaged M. is ap- 
proximately proportional to either ~ s~ 0,1 (H-model) or to 
~ s~ 0,2 (C-model). In the limit of a nearly constant Mach 
number, wave perturbations assume the form of Archimedes' 
spirals. 

As pointed out in Paper II, in a three-dimensional space 
circumplanetary spirals are weakened (see also Bate et al. 
2003). In fact, because of the flow circulation above the disk 
midplane (see Paper II, Figure 3, middle and bottom pan- 
els), less angular momentum is carried inside the Roche lobe 
by the midplane flow. As anticipated above, the spiral density 
features look smoother in the models presented here. This is 



quantitatively demonstrated in Figure 17. In the Figure the 
average specific angular momentum s w rot obtained from the 
M p — 1 Mj C-model (solid line) is compared to that of a 2D 
Jupiter-mass local-isothermal model (short-dash line) and to 
that of a 3D Jupiter-mass model from Paper II (long-dash 
line), which also assumes a quasi-isothermal gas around the 
planet. In the 3D case, the quantity sw TOt is computed at 
z = 0, i.e., in the disk midplane. The latter two models 
are quasi-isothermal in the sense that the temperature dis- 
tribution is fixed and depends only on the distance distance 
from the star r according to the relation T oc (p/7) h 2 M*/r. 
Thus, around the planet's location r = r p , T w 90 K if the 
disk aspect ratio h is 0.05, n = 2.39, and 7 = 1.4. From 
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Fig. 15. — Density distribution and velocity field around the Jupiter-mass model in which the a-viscosity prescription is used 
and set to a = 1CP 2 (left) and a = 1CP 3 (right). The color bar scale is linear and 10 -4 is equal to £ = 32.9 gem -2 . 



Figure 17 one can see that in the outer part of the Roche 
lobe (s > 0.3 -Rh), where spiral perturbations are strongest, 
the discrepancy between the two quasi-isothermal models lies 
between 20 and 30%, whereas the two-dimensional thermal 
model provides a specific angular momentum only 10% larger 
than that of the three-dimensional quasi-isothermal model. 



s/R n 

■J 1 Q.2 05 4. 0.5 D.fi 
I 



■ 



i .DOG rV 



QAM 



0.C1C 



'i 




(J. 10 0.15 
* \AU] 



Fig. 16. — Surface density averaged around M p = 1 Mj plan- 
ets in the inner part of the Roche lobe. Considered models 
are H (solid line), C (short-dash line), A2 (Dash-dot line), 
and A3 (long-dash line). 
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Fig. 17. 

around M ] 

and thermal structures. The solid line belongs to the 2D C- 
model presented in this paper. The short-dash line derives 
from a locally isothermal 2D model (T ~ 90 K) as described 
in Paper I. The long-dash line comes from a locally isother- 
mal 3D model (T ~ 90 K) presented in Paper II. In this last 
case, the specific angular momentum is evaluated in the disk 
midplane (z = 0). 
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Fig. 18. — Toomre parameter Q = c s w ro t/{^ G s S) aver- 
aged around M p = lMj bodies in models H (solid line), C 
(short-dash line), and A3 (long-dash line). The large value of 
(Q) implies that circumplanetary disks around Jupiter-mass 
planets are gravitationally stable. 



This occurrence can be attributed to the better pressure 
structure achieved by the calculations that account for the 
thermal energy budget. 

Circumplanetary disks should not suffer from gravitational 
instabilities since the Toomre parameter Q is orders of mag- 
nitude larger than unity, as proved by Figure 18. This arises 
from the small value of the ratio between the mass of the 
circumplanetary disk and that of the planet. 

Two-dimensional temperature distributions are rather 
symmetric with respect to the planet's position, though they 
are marked by weak spiral perturbations. Figure 19 (top 
panels) shows the azimuthal averaged temperature around 
protoplanet with M p = 1 Mj (left) and M p = 0.1 Mj (right). 
In the Jupiter-mass case, maximum temperatures range from 
roughly 1500 to 1000 K, whereas T reaches the value of 
the ambient medium (< 50 K) toward the border of the Hill 
sphere. Comparing H- and C-models, one can realize that the 
average temperature in the inner part of the Roche lobe differ 
by less than a factor two, whereas the a-viscosity model A3 
shows significant lower temperatures. In the low-mass case 
(M p =0.1 Mj), the peak temperature is around 1000K and 
between w 80 and w 40 K at s = Rh, respectively for both H- 
and C-model. The average temperature differs by less than 
~ 30 K (Fig. 19, upper-right panel). 

Over the entire sub-disk domain, (T) can be well fitted by 
the curve: 

<T)^<T> 3 - (J)*, (38) 

in which the length s is set to 0.1 Rh, for convenience. Fitting 
parameters are reported in Table 4, together with the validity 



range of the fitting function. From entries in the Table, it 
turns out that the temperature generally falls off as 1/s. 

Bottom panels of Figure 19 display the average effective 
temperature (T e a) (see eq. [11]) in the vicinity of protoplan- 
cts. Figure 20 shows that the magnitude of the optical thick- 
ness (t) is between 10 and 100, hence equation (14) yields 
T e ff ~ jy r °- 28 anc j thus there is, at most, a factor three be- 
tween the two temperatures. It should be noted that only the 
model A3 presents optically thin regions toward the edges of 
the sub- disk. 

Accreted matter also contributes to the energy budget of 
the circumplanetary disk via dissipation of gravitational en- 
ergy into heat. This source of energy is not treated explicitly 
in these simulations. However, it is considered implicitly as 
a compression work term arising from the removal of mat- 
ter in the accreting region. Tanigawa & Watanabe (2002) 
assumed the dissipation of gravitational energy as the only 
heating source and derived temperature profiles of sub-disks 
with different sound-speed regimes, performing local simula- 
tions with an isothermal equation of state. They found that 
the temperature would scale as 1/ ^fs, independently of the 
planetary mass. 

We have to point out that, since radial transport of radi- 
ation is not treated here, one should expect large tempera- 
ture gradients to be smoothed out. But we have seen that 
temperature profiles are not extremely steep, in the range of 
distances covered by these models, hence correction by radial 
transfer should not be dramatically relevant, as we are going 
to demonstrate. 
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Fig. 21. — Average circumplanetary disk aspect ratio 
around M p = 1 Mj protoplancts. The solid, short-dash, 
and long-dash lines reproduce the H-, C-, and A3-model, 
respectively. 
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Table 4 

Fit Parameters for the Averaged Temperature Distribution: Accreting Models. 



Mp/Afj 


C-MODELS 


H-MODELS 


(T)s [Kj C 


Range 


(T)s [Kj £ Range 


1.0 
0.5 
0.2 
0.1 
0.06 


5.405 x 10 2 0.97 
5.304 x 10 2 0.96 
5.713 x 10 2 1.12 
4.395 x 10 2 1.03 
4.634 x 10 2 1.07 


[0.05,1.0] Rh 
[0.05,1.0] Rr 
[0.05,1.0] R R 
[0.05,1.0] i? H 
[0.05,1.0] i? H 


9.534 x 10 2 0.96 [0.1, 1.0] Rn 
8.425 x 10 2 0.95 [0.1,1.0]i? H 
6.713 x 10 2 0.90 [0.05, 1.0] Ru 
4.868 x 10 2 0.79 [0.05, 1.0] i? H 




A3-MODEL 


A2-MODEL 


1.0 


2.484 x 10 2 1.04 


[0.05,0.8] Ru 


4.234 x 10^ 0.95 [0.1, 1.0] Ru 



Note. — Parameters that enter equation (38). This is a linear best-fit of log(T), 
computed inside of the Hill sphere. The column "Range" refers to the validity range 
of the fit. 



6.1.1. Radial Radiation Transfer in Circumplanetary 
Disks 

Since the temperature gradient in the vicinity of a proto- 
planet may be quite steep, one of the hypothesis used in § 2.1 
might locally break down. We refer to the assumption that 
the radiative flux in the vertical direction overwhelms those 
in the horizontal direction. With a semi-analytical analysis 
it is possible to make an a posteriori check of that hypoth- 
esis and shed some light over this matter. For simplicity we 
will suppose that the temperature distribution has a cylin- 
drical symmetry around the planet, which is indeed what is 
observed in all of the models. Hence the total flux of energy 
transferred via radiation is 



A ■ 



J-oc dz J_ ; 



ld_ 

s ds 



s F s ] dz = A z + A s 



(39) 



where s is the distance from the planet. The two terms on 
the right-hand side can be written as 



A z ~- 



128 ctr r 
3«p H ' 



and 



6 k p s z 



(40) 



(41) 



Equation (40) was obtained by setting dT/dz « T/H in 
equation (39), whereas equation (41) was derived by adopt- 
ing equation (38), which we have checked to hold as long as 
s < s < Ru- The validity of the relation (9), namely that 
A ~ A z , depends on the ratio of the left-hand sides of equa- 
tions (40) and (41): 



| A s | _ 2 
|A,| ~ ? 



H 



(42) 



According to equation (42), only when H < s the radia- 
tive cooling (eq. [14]) included in the energy equation is also 



locally a good approximation. As shown in Figure 21, this 
seems indeed to be the case in models involving Jupiter-size 
planets. To cast the above ratio in a more explicit form, one 
can use the expression for H (eq. [19]) in the limit s <C r, 
which yields 



H 



£s 
7 



G Mr 



(43) 



By using the form of the sound speed in equation (5), one 
gets 

" V ' kT < (44) 



v /i toh J \G Mp y 

Therefore, because of the assumption made above on the tem- 
perature profile, equation (42) becomes 



\A4 = 2 / fc(T)s 
|A Z | \ firriH 



G Mr: 



(!) 



5-1 



(45) 



Since the distance s is constrained by the applicability range 
of equation (38), the ratio given by equation (45) is suppos- 
edly meaningful roughly for s £ [0.1 -Rh,-Rh], for the models 
considered here. With the appropriate values, the ratio (45) 
becomes 



(T)s 
1000 K 



M, 

M v 



2/3 



5-1 



(46) 



Equation (46) is an evaluation of equation (45) at a distance 
s from the planet, whose length shortens as Mp 1 / 3 . Since 
the squared sub-disk aspect ratio grows as 1/M P , the flux 
ratio goes as M p 2 ^ . From Table 4, we see that £ — 1| is 
between 0.03 and 0.12 in C-models and between 0.04 and 
0.21 in H-models. Therefore, the right-most term in equa- 
tion (46) can be considered as a unity term. Thus, around 
s ~ s = 0.1 Ru, A s |/]Az| = 0.13 for the Jupiter-mass H- 
model and a factor two as small for the C-model. Similar 
small numbers are obtained from both the A2 and A3 model. 
For Mp/Mj — 0.1, the ratio increases, respectively, to 0.2 
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Fig. 19. — Azimuthally averaged temperature (top) and effective temperature (bottom) distribution, around Jupiter-size (left) 
and 0.1 Mj (right) accreting models. Solid, short-dash, and long-dash lines refer to models H, C, and A3, respectively. 



and 0.3, which is consistent with the applicability limit of the 
two-dimensional approximation. Thereby, apart from regions 
closer than ~ 0.1 Rn to the accreting planet which, probably, 
cannot be properly investigated by these computations, radi- 
ation transport in the radial direction does not play a major 
role in the energy budget of circumplanetary disk material. 
This is in agreement with the standard scenario of the late 
stages of protojovian disks (e.g., Coradini et al. 1989; Canup 
& Ward 2002). 

6.2. Non-accreting Models 

Accreting models assume that the planetary accretion rate 
can adjust to the rate at which matter can be supplied by the 
surrounding environment. We now investigate the opposite 



situation, namely when no further accretion onto the proto- 
planet is allowed. 

The overall flow structure and density distribution of non- 
accreting protoplanets are drastically different from those of 
accreting ones. This can be clearly seen in Figure 22, where 
we show the surface density and the flow field for M p = 1 Mj 
(upper panels) and 0.1 Mj (lower panels) objects. The stan- 
dard scenario of a Keplerian-like sub-disk is completely al- 
tered. This is a consequence of the large pressure gradient, 
built up by the large density gradient, which allows matter 
to move on orbits not constrained by the centrifugal balance. 
Therefore, the two-arm spiral feature is replaced by a more 
complex system of multiple shock fronts across which material 
is first deflected toward the planet and then ejected outward, 
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Fig. 20. — Average optical thickness around M p = 1 Mj (left) and M p = 0.1 Mj (right) protoplanets. The solid, short-dash, 
and long-dash lines reproduce the H-, C-, and A3-model, respectively. 



as also indicated by the sign of the in-fall velocity iOi n within 
the region. 

The flow strongly diverges in a zone that has a radius equal 
to « 0.4 Rh in the upper-left panel (H- Jupiter model) and to 
~ 0.2 Rh in the upper-right panel (C- Jupiter model). The 
average surface density inside of Rh/3 is between one and 
two and a half orders of magnitude larger than that mea- 
sured in accreting models. The mass collected in s < 0.8 Rh 
is 1.97 x 10~ 4 Mj and 4.30 x 10~ 5 Mj for the Jupiter-mass H- 
model and the C-model, respectively. The azimuthal average 
of the circumplanetary density, around Jupiter-size objects, 
can be approximated to (S) ~ 1.2 x 10 4 (0.1 Rh/s) 2 ' 8 gem -2 
(H-model) and to (£} ~ 1.5 x 10 3 (0.1 Rh/s) 2 ' 7 gem" 2 (C- 
model). The top panel of Figure 23 shows the averaged sur- 
face density profiles around M p = 1 Mj (left). 

A major difference can be also observed in the sub-disk 
circulation of low-mass models (Fig. 22, lower panels). In 
fact, the fluid rotates in a clockwise direction. In general, 
the direction of rotation in sub-disks is determined by the 
balance among the Coriolis force, the pressure gradient, and 
the gravitational attraction by the planet. Basically, refer- 
ring to the lower-right panel of Figure 22, the Coriolis force 
deflects fluid elements, orbiting at r < r p , rightward. There- 
fore, it forces matter to cross the gap and to reconnect to the 
other side, while still moving in a position upstream of the 
perturber body. In contrast, the term — dP/ds, supposedly 
positive, opposes to reconnection when fluid is upstream of 
the planet's location but favors it when matter is downstream 
of the planet. The prograde rotation that has been encoun- 
tered so far indicates that the Coriolis deflection overwhelms 
the pressure gradient. For this reason reconnection across 
the gap always occurs upstream of the planet. Evidently, in 
the low- mass models (M p = 33 Mg) displayed in the bottom 
panels of Figure 22, the pressure gradient drives material, 
flowing from upstream, past the planet's position and forces 



it to circulate clockwise around the perturber. Consequently, 
matter entering the Hill sphere has an angular momentum 
anti-parallel (in the planet's frame) to that of the circumstel- 
lar disk. In the range from 0.1 to IRh, the magnitude of 
the rotational velocity |(w ro t)| increases roughly as yfs. The 
behavior of (S) resembles that of a power-law, as in Jupiter- 
mass cases, with powers 0.9 and 0.5, respectively for the H- 
and C-model (Fig. 23, right panel). Interestingly, no large 
differences are seen in the two density profiles. 

The high pressure gradient is also caused by the large tem- 
perature gradient. Parameters for analytic approximations of 
the average temperature (T), according to equation (38), are 
reported in Table 5. Maximum temperatures, at 

s = 0.1 Rh, reach 1900 K and 1200 K, in M p = 1 Mj H- and 
C- models respectively (bottom-left panel in Fig. 23). At the 
limit of the Roche lobe the temperature is somewhat higher 
than that measured in accreting models and depends on the 
gap structure. They vary from ~ 100 to ~ 200 K. Low mass 
models furnish temperatures that are very similar in both 
viscosity regimes. 

It is worthy to notice that in either planetary-mass and 
viscosity case, the resulting temperature profiles are steeper 
than the ones shown in Figure 19 only when s < 0.1 Rh, but 
they are nearly as steep past this distance. Actually, from 
Table 5 it appears the toward low planetary masses temper- 
ature profiles are even flatter. Performing the same check as 
done in the previous section, it turns out that |A S | / |A Z is at 
most 0.3. In fact, around Jupiter-mass planets, the ratio H/s 
varies between 0.3 and 0.5, for distances larger than 0.1 _Rh- 
Hence, even in these circumstances the type of energy equa- 
tion implemented here yields a reasonable description of the 
thermal structure of sub-disks, down to ~ 0.1 Rh- 
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Fig. 22. — Density distribution and velocity field around M p = 1 Mj (top) and M p = 0.1 Mj (bottom) non-accreting models. 
Left panels show H-models while right ones show to C-models. Colors scale logarithmically in all of the panels. In the units 
used in the color bar, —3 corresponds to £ = 328.7 gcm~ 2 . 



7. Accretion and Migration 

The evaluation of the mass accretion rate of protoplanets is 
carried out by reducing the surface density around the planet, 
according to a relation of the type AE/S = At/r cv (see, 
Paper I for details). This reduction process has a time scale 
r cv that is much smaller than the integration time step At. 
In Figure 24 the planetary accretion rate M p evaluated in 
these computations is compared to that achieved with local- 
isothermal 2D (Paper I) and 3D (Paper II) models, where 
the disk mass was rescaled to match the values adopted here 
(see § 4). 

In the context of the accretion history of protoplanets, it 



is worthwhile to stress that 2D calculations are only able to 
simulate the late stages of the accretion process, i.e., when it 
mainly occurs via mass transfer in the circumplanetary disk. 
This is consistent with the parameterization adopted here to 
describe the circumstellar environment. 

From Figure 24, one can see that C-models (filled squares) 
provide accretion rates quite similar to those obtained from 
models with a fixed temperature distribution (crosses). This 
occurrence is related to the quasi-Keplerian, circumplanetary 
flow observed in both kinds of computations. Jupiter-mass 
models with an a-type viscosity (open triangles) give rates 
which are similar to that of the C-model when a — 1CP 2 but 
somewhat lower (M p = 6 x 1CT 3 M e yr _1 ) when a = 1CT 3 , 
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Fig. 23. — Average surface density (top) and temperature (bottom) around 1 Mj (left) and 0.1 Mj (right) non-accreting 
protoplanets. Solid lines refer to the H-models whereas short-dash and long-dash lines belong to C-models. 



consistently with the low sub-disk mass measured in this case. 
H-models provide values of M p (filled diamonds) that, in case 
of Jupiter-size planets, fall above the estimates from all of 
the other models. According to the accretion procedure, the 
larger the amount of matter contained in the accretion region 
s = K ac , the larger the accreted mass is. Within s — 0.1 Rn of 
a Jupiter-mass planet, the mean density computed in the H- 
model is roughly three times larger than that in the C-model 
(Fig. 14). This is related to the background gap density that 
is much higher in H-models, as reported in Table 2 (see also 
Fig. 7), due to the larger fluid viscosity. 

Figure 25 shows the total gravitational torque 7d exerted 
by disk material onto accreting planets, as the inverse of the 
migration time scale (7d oc \r p \/r p = 1/tm). This quantity is 



computed from the gravitational force acting on the planet, 
as explained in Paper I. The length j3 represents the ra- 
dius of the circular region around the planet's position whose 
contribution is not taken into account. In Figure 25, this 
length is set equal to 7?h- Outcomes from C-modcls, and 
generally those from H-models, compare well to the local- 
isothermal computations. Models with a-viscosity provide 
migration rates on the same magnitude as well. C-models 
show a rather constant value of tm as a function of the planet's 
mass, whereas H-models present a minimum at M p = 0.2 Mj. 
Migration time scales of H-models differ by less than a factor 
two from those of C-models. Only for M p = 33 M®, tm is 
appreciably longer. One may attribute this effect to the ab- 
sence of a real gap in this model (see Table 2), and thus to 
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Table 5 

Fit Parameters for the Averaged Temperature Distribution: Non- Accreting Models. 



Mp/Afj 


C-MODELS 


H-MODELS 






Range 


CO,- [k] 




Range 


1.0 


1.187 x 10 3 


1.29 


[0.05,1.0] Rn 


1.909 x 10 3 


1.11 


[0.05,1.0] i? H 


0.5 


1.132 x 10 3 


0.81 


[0.05,1.0] i? H 


1.458 x 10 3 


1.00 


[0.05, 1.0] i? H 


0.2 


9.256 x 10 2 


1.18 


[0.05, 1.0] i? H 


9.335 x 10 2 


0.76 


[0.05, 1.0] i? H 


0.1 


7.570 x 10 2 


0.70 


[0.05, 1.0] i? H 


6.242 x 10 2 


0.59 


[0.05, 1.0] i? H 


0.06 


6.103 x 10 2 


0.60 


[0.05, 1.0] i? H 









Note. — Parameters that enter equation (38) when applied to non-accreting models. 
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Fig. 24. — Comparison of the accretion rates obtained from 
different types of simulations. Filled diamonds and squares 
indicate estimates provided by H- and C-models. Open tri- 
angles show data points from a-viscosity models A2 (larger 
value) and A3. Asterisks and crosses represent the rates fur- 
nished by local isothermal computations in two (Paper I) 
and three (Paper II) dimensions, respectively. These last 
two sets of values are rescaled to match the disk mass Md = 
4.8 x 10 -3 M«, used in these simulations. 



gap region is fixed to 5 Rn (see Fig. 7) and the zone behind 
the planet (&) is distinguished from the one ahead of it (C). 
Torques exerted from inside the Hill sphere are not accounted 
for. A large fraction of the contributions arising from zones 
B and C is actually generated in the coorbital regions, within 
a few Hill radii from the planet. From the entries in Table 6 
one can argue that coorbital torques in H-models are larger in 
magnitude than they are in C-models. This is a consequence 
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an early onset of the Type I drifting regime. 

[!t] In Table 6, we separate the contributions of different 
disk regions to the total torque 7r>, introducing the inner 
disk (A), the outer disk (Z>) and the gap region (£> and C). 
Further details on the regions are reported in the Table's note. 
All contributions are normalized to |7d|. The width of the 



Fig. 25. — Migration time scales plotted against the plan- 
etary mass. Notations are the same as in Figure 24. Data 
points indicate the magnitude of the total gravitational torque 
when this is computed excluding the contributions of matter 
residing inside of the sphere s — Ru ■ Only accreting models 
are considered. As for a-viscosity models (open triangles), 
A2 gives the faster migration rate of the two. 
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Table 6 

Torques From Different Disk Regions: Accreting Models. 



Normalized Torques 

Model 



A B C V 



1.0 Mj 


c 





.0987 


-1.9416 


1.3707 


-0. 


.5278 


1.0 Mj 


H 





.3247 


-13.538 


12.828 


-0 


.6147 


1.0 Mj 


A3 





.9857 


-4.7780 


4.9240 


-2. 


.1317 


1.0 Mj 


A2 





.0913 


-2.1357 


1.3905 


-0. 


.3461 


0.5 Mj 


C 





.1962 


-3.6112 


2.9099 


-0. 


.4949 


0.5 Mj 


H 





.3301 


-14.940 


14.084 


-0 


.4732 


0.2 Mj 


C 





.4280 


-13.852 


13.176 


-0 


.7520 


0.2 Mj 


H 





.2406 


-14.827 


13.840 


-0. 


.2536 


0.1 Mj 


C 





.4070 


-17.989 


17.211 


-0. 


.6290 


0.1 Mj 


H 





.6315 


-36.823 


35.768 


-0. 


.5765 


20 M e 


C 





.3471 


-17.656 


16.800 


-0. 


.4911 



NOTE. — Gravitational torques exerted by the inner 
disk (A) r < r p — 2.5 Rn, and the outer disk (V) r > 
r p + 2.5 Rh- The remaining domain (r p — 2.5 Ru < r < 
r p + 2.5 Ru) is further divided in the trailing gap region 
(B) ip < ip p , and the leading gap region (C) ip > (p p . The 
width of the gap region is chosen according to the average 
width observed in Figure 7. All torques are normalized to 
the absolute value of the total torque T D . Furthermore, 
they exclude the contribution from a circle of radius Rn 
centered on the planet. 
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of the shallower gap of H-models. The situation is reversed 
in the a-viscosity models because of the density humps at L4 
and L5 Lagrangian points. 

As for the influence of the parameter f3 on the total torque, 
it turns out that H-models are generally more sensitive than 
C-models. By examining the function 7d = To(/3), some 
more insight can be gained. For M p = 1 Mj models (Fig. 26, 
top- left panel), no significant net torque arises from within 
the Hill sphere in the C-model (dashed line). The situation 
appears to be different for the H-model (solid line) because 
Td decreases to a minimum, around f3 = Rh/2, and then 
increases and changes sign around /3 = 0.1 Rh- Jupiter- 
mass a-viscosity models (Fig. 26, top-right panel) display 
a smooth behavior since the total torque is rather constant 
down to /3 = -Rh/2. When (3 diminishes, To increases (i.e., 
it has smaller absolute values) and the model with a — 1CF 3 
(dashed line) shows a change of sign at /3 = 0.4 Rh- When 
M p =0.1 Mj (Fig. 26, bottom- left panel), in either case, 
dTo/d/3 < 0, because of the positive torques arising within 
the sub-disk. The rate of change is more drastic in the H- 
model, where the total torque becomes positive and keeps 
increasing when j3 < -Rh/2. Torques acting on the 20 
planet (Fig. 26, bottom-right panel) are strongly dependent 
on f3, as positive torques are exerted down to w 0.4 Rh and 
large negative torques arise between « 0.2 Rh and m 0.4 Rh. 

8. Conclusions 

In this paper we presented two-dimensional simulations 
of disk-planet interactions in which we considered the joint 
thermo-hydrodynamics evolution of the system. In order to 
do that, we solved an energy equation that accounts for the 
major processes responsible for the energetic balance of fluid 
parcels. We restricted the treatment of the problem to a 
flat geometry in order to make some simplifying assumptions 
on the radiative part of the energy budget, i.e., to treat the 
radiation transfer as a cooling term. A nested-grid technique 
was used to resolve different length scales, therefore both the 
global and local features could be investigated. 

The two-dimensional approximation has been proved to 
work reasonably well as long as the planetary mass is larger 
than a tenth of the mass of Jupiter (Paper II). Because of 
this, only planets more massive than « 0.1 M p have been 
simulated. In order to compromise with the lingering uncer- 
tainty on the viscosity magnitude and to inquire into different 
temperature regimes, we considered different constant viscos- 
ity and a-viscosity models. Both accreting and non-accreting 
planets were simulated. 

From the global point of view, and roughly speaking, we 
can conclude that circumstellar disk models with a fixed tem- 
perature distribution provide a reasonably good description 
of the system, though details of the gap structure and shape 
of the disk spirals differ. Around Jupiter-mass protoplanets, 
only with kinematic viscosities as low as 10 15 cm 2 s _1 , a wide 
and deep gap is carved in whereas it reduces to a shallow 
density trough when the viscosity is ten times as large. Mean 
temperatures in the gap can be as low as 20-50 K, depending 
on the viscosity regime. The most diluted and cold gap ma- 
terial is located where the circumplanetary disk merges into 
the gap. The aspect ratio appears smaller than 0.04 and it 



is also marked by a gap-like feature, which shades the pro- 
toplanet's environment against stellar radiation. The models 
with a = 10~ 3 show a ring of matter within the gap that 
evolves on a viscous time scale. Regarding smaller planetary 
masses (M p « 0.1 Mj), both gap and spiral perturbations 
become quite faint features. Likely, only Jupiter-size bodies 
are capable of digging gaps deep enough to be observed by 
present-day instrumentation. 

From the local point of view, we obtained the distribution 
of density, temperature, and other thermodynamical quan- 
tities in circumplanetary disks. Except for the a-viscosity 
model with a = 10~ 3 , rather comparable outcomes are pro- 
vided by models with different viscosity regimes. The den- 
sity, averaged around the planet, shows an exponential ra- 
dial decline. Sub-disks have masses on the order of 10~ 5 - 
10~ 6 Mj, which makes them gravitationally stable. Tem- 
peratures range from several hundred Kelvin degrees, at dis- 
tances s < O.IRh from the planet, to values between 50 and 
100 K, at the border of the Hill sphere. Azimuthal aver- 
ages yield profiles dropping nearly as the inverse of the dis- 
tance from the planet. Circumplanetary disks have aspect 
ratios H/s around a few tenths and they tend to be optically 
thick. All of these properties appear to be consistent with the 
"fast-inflow" analytical models of protojovian disks recently 
developed by Canup & Ward (2002). By means of a semi- 
analytical treatment, we showed that, also in the sub-disk, 
the neglect of radial radiation transfer should not crucially 
affect our results. Non-accreting models furnish a rather dif- 
ferent scenario. Around Jupiter-mass objects, circulation in 
the Roche lobe is not Keplerian-like because of the large pres- 
sure built up mainly by the density gradient and partly by the 
temperature gradient. Around M p = 0.1 Mj, clockwise ro- 
tation establishes as a result of the weakening of the Coriolis 
force with respect to the pressure gradient term. 

Accretion and migration rates inferred by these simula- 
tions are comparable to those evaluated with previous purely 
hydrodynamical local-isothermal computations. An analysis 
of the zonal torques indicates that coorbital torques increase 
as viscosity increases. Moreover, negative coorbital torques, 
arising from material outside of the Roche lobe and lagging 
behind the planet, are generally larger than those from ma- 
terial ahead of the planet. We also found that, in the two- 
dimensional approximation, the material lying inside of the 
Roche lobe is able to exert strong gravitational torques on 
the planet. Since the density gap is usually filled in H-models 
(see Table 2), Type I migration regime might extend to larger 
planetary masses (see Fig. 25). 

In order to address the issue of how the third dimension af- 
fects the thermo-hydrodynamics of the system, both globally 
and locally around protoplanets, we are extending this study 
to 3D simulations in which radiation transfer is treated in the 
flux-limited diffusion approximation (Levermore & Pomran- 
ing 1981). 

We thank Udo Ziegler for having made available to us 
an early FORTRAN Version of his code Nirvana. G. D. is 
deeply grateful to all the friends at the INAF-OAC (Naples) 
and RSMAS/MPO (University of Miami) for their warm hos- 
pitality. The paper benefited from the useful comments of an 




Fig. 26. — Inverse of the migration time scale versus the radius (5 of the region excluded from the torque calculation. The 
quantity v = a/\d\ indicates the direction of the planet migration. Solid lines refer to models H or A2 while dashed lines refer 
models C or A3. Top-left panel. Jupiter-mass models H and C. Top-right. Jupiter-mass models with a-viscosity A2 and A3. 
Bottom-left. M p = 0.1 Mj model H and C. Bottom-right. M p = 20 M e C-model. 
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